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Asymmetric nuclear matter is investigated in the low density region below the nuclear saturation 
' density. Microscopic calculations based on the Dirac Brueckner Hartree-Fock (DBHF) approach 

' with realistic nucleon-nucleon potentials are used to adjust a low density functional. This functional 

, is constructed on a density expansion of the relativistic mean field theory which allows a clear 

interpretation of the role of the mesons to the equation of state. It is shown that a correction 
term should be added to the functional in order to take into account the effects beyond the mean 
field. Two functionals with different corrections are obtained and their topological properties has 
been studied. Those functionals converge to predict a reduction of the spinodal zone in asymmetric 
' nuclear matter by about 15-20% and an isoscalar unstable mode closer to the constant Z/A direction 

04 ^ than the functional without correction. 
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How does nuclear matter properties change when the density decreases from saturation densities ? In this low 
density regime, what is the role of the isospin asymmetry ? Indeed, the nuclear density functional below saturation 
density has not been much studied while it has an importance for several topics concerning atomic nuclei surface 
04 ' properties, the equation of state of the core of neutron stars and for the dynamical description of heavy ion collision, 
^ \ both at intermediate and relativistic energies. Recently, several attempts have tried to establish a relation between 
1^ the low density equation of state and the nuclear properties like surface behavior and pairing properties [J, H, d, 0] j 
neutron radii 5] or the spinodal instability 3 . Those works are based on phenomenological density functionals or fits 
of ab-initio calculations but without considering the very low properties of the equation of state. Below saturation 
density, effects of two-body correlations are important and induce anomalies in the density dependence of the equation 
of state ■ 

Models which make predictions for the nuclear equation of state (EOS) can roughly be divided into three classes: 
phenomenological density functionals, effective field theory approaches, and ab initio approaches. The phenomeno- 
i;^ logical density functionals are based on effective density dependent interactions such as Gogny or Skyrme forces [1, 

or relativistic mean field (RMF) models Parameters are adjusted to nuclear bulk properties and finite nuclei. 

Effective field theory (EFT) approaches are based on a perturbative expansion of the nucleon-nucleon interaction or 
^ ' the nuclear mean field within power counting schemes. These approaches lead to a more systematic expansion of 
the EOS in powers of density, respectively the Fermi momentum kp- The EFT approaches can be based on density 
functional theory [lll[l2| or e.g. on chiral perturbation theory The advantage of EFT is the small number 

of free parameters and a correspondingly higher predictive power. However, when high precision fits are intended, 
the EFT functionals are based on approximately the same number of model parameters as phenomenological density 
functionals due to fine tuning through additional parameters. 

Ab initio approaches are based on high precision free space nucleon-nucleon interactions. In addition, predictions 
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for the nuclear EOS a re p arameter free. Examples of such approaches are variational calculations ISLlK 



Brueckner- 



i re p arameter ire' 

Hartree-Fock (BHF) IH, El H or relativistic Dirac-Brueckner-Hartree-Fock (DBHF) [2l|, |2j, l23|, l2j, l25|, |2e 
calculations and Green function Monte-Carlo approaches (GFMC) [27|, [H, [2^ . Non-relativistic ab initio calculations 
do not meet the empirical region of saturation, whereas relativistic calculations do a better job. This deficiency can 
be solved by the explicit inclusion of three-body forces where the relativistic approach accounts already effectively for 
part of these contributions. For a more detailed discussion see e.g. Ref. [20j. 

In the present work microscopic calculations based on the DBHF approach [26| using a realistic nucleon-nucleon 
potential, i.e. the Bonn A interaction [30| . are used to obtain a functional which describes the equation of state from 
low densities u p t o saturation density. The construction of this functional is motivated by relativistic mean field 



rp to sa 

(RMF) theory |3l|, IS^. Since both, DBHF and RMF are relativistic approaches based Dirac phenomenology they 
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have similar features, in particular the occurrence of large and cancellating scalar and vector fields in the isoscalar 
sector. That the occurrence of such fields is a fundamental consequence of the elementary nuclear force has recently 
been shown in [s^l- Similar fields, which are, however, smaller, occur also in the isovector sector. Thus, RMF theory 
is well suited for the present investigations. However, RMF theory is insufficient to reproduce the more complex 
nonlinear behavior of the DBHF energy density near ps =0.1 fm^'^, where for instance effects of the deuteron pole 
show up. Corrections beyond mean field are necessary and the non-linear behavior is then corrected in adding new 
terms in the functional. 

The paper is organized as follows: the relativistic DBHF is shortly sketched in Sec. HI Furthermore, Sec.[n]is devoted 
to the relativistic mean-field approach and the series expansion in the baryonic density. This expansion allows a clear 
interpretation of the meson contributions to the equation of state. The parameters of the functionals are obtained 
in Sec. IIIII With these density functionals the dynamics of the liquid-gas transition induced by heavy ion collision 
at Fermi energies, i.e. the spinodal instabilities, are investigated and analyzed in Sec. IIVI Finally, we end with a 
conclusion in Sec. |Vl 

I. DBHF APPROACH 

We consider homogeneous nuclear matter at low density so that two-body correlations dominate. Of course, at very 
low density, clustering phenomena can occur, like deuteron, tritium, helium and alpha particle formation 0, Is^ . l35l |. 
Therefore, we are going to consider densities which are low with respect to saturation density, but still large compared 
to typical values where the onset of clustering occurs. Typically, we are considering densities between about one 
tenth to one half of the saturation density. In this density region, microscopic calculations based on the relativistic 
Dirac-Brueckner-Hartree-Fock (DBHF) approach [2l|, [23, [2^ [2^, [Sg] are expected to be quite accurate, and they will 
be, therefore, the starting point of our analysis. 
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FIG. 1: (Color online) Comparison of several equation of state with the DBHF results. The DBHF calculation are represented 
by the stars, the variational by the empty squares, the non-relativistic Brueckner results by the dashed line and two different 
relativistic mean field parametrisations: NL3 (dashed-dotted line) and DD-TW (dotted line). 

In the relativistic Brueckner approach the nucleon inside the medium is dressed by the self-energy E. The in- 
medium T-matrix which is obtained from the relativistic Bethe-Salpeter (BS) equation plays the role of an effective 
two-body interaction which contains all short-range and many-body correlations of the ladder approximation. Solving 
the BS-equation the Pauli principle is respected and intermediate scattering states are projected out of the Fermi 
sea. The summation of the T-matrix over the occupied states inside the Fermi sea yields finally the self-energy in 
Hartree-Fock approximation. This coupled set of equations states a self-consistency problem which has to be solved 
by iteration. Technical details of the present DBHF calculations, in particular the treatment of isospin asymmetry, 
can be found in [26| . The results are based on the Bonn A one-boson-exchange potential for the bare nucleon- nucleon 
interaction [30| . 

Fig. ID) compares the prediction for the low density regime of symmetric nuclear matter (below zero) and pure 
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neutron matter (above zero) from microscopic many-body approaches. It shows in addition the equations-of-state ob- 
tained by two typical phenomenological relativistic mean field models. From the microscopic side these are the DBHF 
results [2^ which will be further analyzed, non-relativistic Brueckner (BHF) calculations from [l[ and variational cal- 
culations from [1^. The variational and BHF calculations are based on the AVig Argonne potential. The variational 
calculations include in addition phenomenological 3-body-forces (Urbana IX) and relativistic boost corrections which 
both do, however, iiot play an important role in the low density regime. The phenomenological models are the well 
established NL3 fs?! relativistic mean field (RMF) parametrisation and the RMF model of ^ (DD-TW). The latter 
is a phenomenological version of density dependent RMF theory using density dependent meson-nucleon couplings 
[sil which allows for a larger freedom in the adjustment of the EOS. Both approaches fit finite nuclei properties with 
high accuracy. 

The first what becomes evident from Fig. ([T|) is the remarkable agreement of the microscopic approaches (DBHF, 
BHF, variational) concerning the pure neutron matter EOS. This indicates that both, the interaction and the many- 
body schemes are well controlled in the nn sector at low densities. The reason is on the one handside the large nn 
scattering length and on the other side the lack of clustering phenomena {d, a etc) which make the treatment of 
neutron matter at subnuclear densities less model dependent. In this context it is worth noticing that the microscopic 
calculations (BHF/DBHF, variational) are consistent with the low density limit of 'exact' Quantum-Monte-Carlo 
calculations for neutron matter [27| and with the result of the renormalization group Viowk approach [ioj . 

The situation seems to be different for symmetric nuclear matter. The Brueckner calculations show significantly 
more binding than the variational calculations of [16]. On the other hand, the DBHF and BHF results are very close 
and exhibit the same low density behavior: in contrast to RMF theory and also to [l^ one can observe a non-linear 
convergence to zero when the density decreases. This fact is associated to the deuteron channel and can possibly be 
interpreted as a manifestation of the onset of the super-fluid phase. In this context it is interesting to note that a 
recent study of low density nuclear matter Q, based on a virial expansion which includes protons, neutrons and a- 
particle degrees of freedom, revealed a low density EOS which is in qualitative agreement with the DBHF predictions. 
In the virial approach the binding energy goes smoothly to zero for neutron matter while the energy per particle E/A 
minus the (free) kinetic energy in symmetric matter remains practically constant at a value around -8 MeV down to 
extremely low densities (ps — 0.0002 fm^'^) before it rapidly drops to zero (Fig. 15 in 7]). Subtracting from the 
DBHF result the kinetic energy of a non-relativistic Fermi gas 3fc^/10M yields at kp — 0.5 fm^^ {pB = 0.0084 fm^"^) 
a values of -8.4 MeV which coincides remarkably well with the virial low density limit. 

II. SERIES EXPANSION OF THE RMF LAGRANGIAN IN THE BARYONIC DENSITY 

A Lagrangian density of interacting many-particle system consisting of nucleons, isoscalar (scalar cr, vector to), and 
isovector (scalar 5, vector p) mesons is the starting point of the relativistic mean field (RMF) theory. 



where cr is the cr-meson field, tj^ is the w-meson field, is p meson field, (5 is the isovector scalar field of the (S-meson. 

F^y = df^ujy — diyUJf^, G^u = d^pi, — di^Pfj,, and the U{a) is a nonlinear potential of a meson : U{a) = ^o,'^^ + jba"^. 
Dynamical equations deduced at the mean field approximation are presented in the appendix |^ We refer to the 
appendix for all the standard definitions. Hereafter, we introduce the coupling constants fi = gp/rrii for i=(T, S, u and 
p, and the non-linear constant /^*' — a{fa-/ma-)^. In the following the energy density e is expanded up to the power 4 
in proton and neutron densities. The expression of the density of energy is given in the appendix [A] For the linear 
version of the RMF model, i.e. without a nonlinear cr-meson potential U{a) such an expansion can be found in [4]| . 
Here we extend this expansion to the non-linear case and to the isospin sector, i.e. to pp ^ pn (the isovector mesons 
p and 8 are included). Notice however, that in the present form, only the non- linear term with the coupling constant 
a is included because the one with the coupling constant h contributes to higher terms in the density expansion. 
The scalar field is the solution of the following self-consistent equation 

9aO = HPs ""—{gaOf - ^—^{gaof . (2) 
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FIG. 2: (Color online) Representation of the different terms in the serie expansion of the energy functional with respect to 
the baryonic density, up to 0.3 fm~^. The first 3 subsets show the fast convergence in symmetric nuclear matter (SNM) and 
the last subset (d) show the contribution of the isovector mesons in pure neutron matter (PNM). To give an idea of the fast 
convergence, we have multiplied some of the terms by huge constants. The coupling constants (set A NLp5) obtained by Liu 
et al. [421 has been used. 



A low density approximate solution is presented in the appendix |B] The solution is expressed as a function of the 
scalar density ps, 



g^a ^ J^ps - Z pi + o{p%) . 
The low-density expansion of the scalar density p^i is then required. It yields 
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where 7 is the degeneracy of the system. Note that this expansion is also a relativistic expansion in the parameter 

kF^|Ml. 

Neutron and proton Dirac masses, also called the scalar masses, are expressed in terms of the scalar and isoscalar 
fields. Using Eq. ([3]) and Eq. (jl]). The low density expansion of the Dirac masses is given by (- proton, + neutron). 
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Now, we evaluate the full density functional. 
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where ekin is the kinetic term. The density functional is decomposed into several terms: 

f^{Pn,Pp) = MpB + epGiPruPp) + ^LiPn, Pp) + ^NL{Pn,Pp) + o(p*) 



(5) 

(6) 
(7) 
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where the term e^G is the contribution of the free Fermi gas without the rest mass, the term cl is generated by 
the interactions and the Dirac mass and the term e^vL is the correction coming from the non-hnear a couphng. 
In the following we give explicitly the form of those terms, classified according to the power in Fermi momentum 
(=power in density divided by 3) in order to have integer index. The pure kinetic contributions (FG) is epc — 
£fg,5 + eFG,7 + £fg,9 + CFG, 11 up to the power 4 in the densities where 
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The contribution of the mesons (with only linear couplings) and Dirac mass contribution is — eL_6+eL,8+ei,io+eL,ii 
where 
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Finally, the first order corrections induced by the non-linear a-coupling is cnl — ^nl,9 + where 
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The convergence of this series expansion is checked in symmetric nuclear matter (SNM) and pure neutron matter 
(PNM) using the set of coupling constants, set A NLpS, obtained by Liu et al. [13] . We choose this set of parameters 
because it has been obtained with the same degrees of freedom as the one we consider in our Lagrangian. We show 
in Fig. [5] the contribution of the different terms of the functional up to the power 4 in the densities. Some terms has 
been multiplied by a huge factor to distinguish from each others. The convergence is essentially due to the shorting in 
power of kp/M which comes with our expansion. This figure shows how negligeable are the terms with large power 
counting in the density expansion, even for the higher densities represented (about 0.3 fm~'^). 

In the following the low density RMF functional defined in Eqs. ([7))- (fT71) will be used to fit the result of the DBHF 
calculation [26| and to calculate the spinodal instabilities in low density asymmetric nuclear matter. 



III. DETERMINATION OF THE PARAMETERS OF THE FUNCTIONAL 



In this section, we fit the result of the DBHF calculation (Dirac mass, energy density and binding energies) in 
the density region between 0.01 and 0.2 fm"'^ using the low density RMF functional defined in Eqs. ((7|- (fT7|) . The 
fits are based on the density expansion of the RMF Lagrangian IT} which contains non-linear terms in the scalar 
a field and linear terms in the vector field oj as well as in the isovector p and 5 fields. An alternative would be 
to perform such a p arametrisation of the DBHF results in terms of density dependent relativistic hadron theory 
(DDRH) [il, [S^ [43, S, nil where non-linearity's due to higher order density corrections are absorbed into density 
dependent meson-nucleon vertices at the level of the effective Lagrangian. The reason why we have chosen the standard 
RMF Lagrangian ([1]) is twofold: this allows a well defined low-density expansion while the density dependence of 
effective meson-nucleon vertices in DDRH depends on the choice of an particular ansatz for these functionals. Secondly, 
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the extraction of such couphng functions from the present DBHF self-energies pG'j shows that the DBHF vector self- 
energy, except for the very low density regime, has a linear density dependence which can be expressed by a linear lo 
meson field. Non-linearities in the scalar channel are absorbed in the non-linear <j terms. The isovector dependence 
can also reasonably well be fitted through the two isovector p and b mesons. In summary, such a procedure allows 
a well defined comparison of the microscopic DBHF model to RMF phenomenology and a controlled investigation of 
the low density regime, where the RMF fits break down and require additional correction terms, as will be seen in 
the following. 

The adjusting procedure is twofold: firstly, we fit the parameters of the RMF Lagrangian using the relativistic Dirac 
mass and the energy density in symmetric and asymmetric nuclear matter obtained from the DBHF calculation. 
The fit includes 23 calculated points between yOB=0-02 and 0.13 fm~'^, plus two densities, pb=0.1658 fm~^ and 
Pb=0.197 fm~'^, for y = pp/ pB=Q to 0.5 with a step=0.05. We obtain the set of parameters called RMF presented in 
Tab.m by fitting the Dirac mass and the density of energy and imposing that the functional passes exactly through the 
point at pb — 0.197 fm~'^ in symmetric nuclear matter and pure neutron matter. The latter condition is imposed in 
order to get a value of the symmetry energy close to DBHF, as shown in Tab.Ull Then, we extract and fit the residual 
difference between the DBHF calculation and the energy per particle in symmetric nuclear matter and pure neutron 
matter separately. In order to check the sensitivity of the results on the functional correction, we have investigated 
two different functionals. In the following, we give the details of the adjusting procedure. 



fits 


U 

[MeV-i] 


[MeV-^W] 


h 

[MeV~^] 


[MeV-^] 


/p 

[MeV^^] 


RMF 
set A NLp5 [42] 


1.693 10"^ 
1.629 10"^ 


3.735 10"* 
9.35 10-* 


7.242 10"^ 
8.013 10"=* 


1.299 10"^ 
1.18 10"^ 


8.843 10~^ 
8.996 10"^ 



TABLE L The parameters of tlie functional reproducing the scalar mass and the energy density of the DBHF calculation is 
compared to the parameters proposed in Ref [42| . 



A. Determination of the a and 5 coupling constants 

The low density expansion of the Dirac mass, Eq. ([5]), is used to determine the linear sigma coupling constant /a-, the 
scalar iso- vector (5-meson fs and the non- linear sigma coupling constant /^'' . We deduce the value of the parameters 
for the adjustment to DBHF results in asymmetric nuclear matter with y = pp/ pB=Q-5, 0.3 and 0.0. In Fig. [3] we 
compare the DBHF results with the best fit. The linear contribution comes from the term f^pB in Eq. then comes 
the quadratic term f^'p%, the isospin asymmetry is essentially coming from the first term /fps the contributions of 
the other terms are negligible. The parameters are given in Tab. U and are compared to the set A NLpS proposed in 
Ref [4^ 1 . The parameters fa and fg are very similar for the two set of parameters while the parameters differ by a 
factor 3. In our fit, /^'' is obtained from the quadratic density dependence of the scalar mass while in the set A NLp6, 
the non-linear a coupling is adjusted to reduce the compressibility modulus. As a consequence, we obtain a lower 
value for the parameter /^', and the compressibility modulus is larger than expected. Moreover, to obtain a good fit 
of the binding energy, the non-linear term eNL,9 defined in Eq. (|16p has to be divided by a factor 2. This may indicate 
that higher order non-linearities in the a field and probably also non-linear lo terms should be taken into account to 
obtain a proper description of the DBHF equation of state beyond saturation density. However, to keep the formalism 
as simple as possible we stick to the strandard NL model and apply this phenomenological correction. The saturation 
properties are shown in Tab. [Hi This indicate that one cannot reproduce the scalar mass density dependence and the 
saturation density with standard a-non linear RMF Lagrangian. This illustrate the convenience of using a density 
functional where the effects of the meson couplings are tractable. In our Lagrangian, the compressibility modulus will 
be lowered by the correction terms induced by the physics beyond the mean field. 

B. Determination of the lu and p coupling constants 

In contrast to DBHF theory, which shows a non-linear convergence to zero in the binding energy, RMF theory 
converges smoothly. We show in Fig. 0] (top panel) a comparison between the DBHF calculation and the low density 
functional RMF for the density of energy, E/V, and the binding energy, E/A. We remind that the non-linear term 
£nl,9 defined in Eq. ^6]) has to be divided by a factor 2. The low density effects are reduced on the energy density 
plot compared to the effects on the binding energy. In fact, as E/V = pbE/A, the low density behavior of the binding 
energy is strongly reduced by the factor ps in the energy density. Therefore, we adjust the mean field functional 
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FIG. 3: (Color online) Comparison between the best RMF adjustment and the scalar mass calculated with DBHF approach. 
The asymmetry parameter y = pp/pB, ranges from 0.5 to 0. It clearly shows a linear (faps) and a quadratic (/^'Ps) behavior 
in the baryonic density. The isospin asymmetry is also well reproduced by the linear term fsps in the asymmetry density ps. 




FIG. 4: (Color online) Comparison of the DBHF results (circles points) and the low density functionals in asymmetric nuclear 
matter. On the top panels are shown the energy density and the binding energy for the parameters RMF, while on the bottom 
panels are shown the binding energies for the corrected functionals RMF+Cl and RMF+C2. 

on the energy density obtained from the DBHF calculation in asymmetric nuclear matter, where the low density 
effects are weaker. We also force the functional to reproduce exactly the DBHF calculations at pB = 0.197 fm~^ in 
symmetric nuclear matter and pure neutron matter. This constraint ensure that the symmetry energy at saturation 
density is well reproduced. The table |I] resumes the parameters obtained from the low density functional RMF. In 
the following section, we propose a correction of the functional to take into account the correlations beyond the mean 
field. 



C. Corrections to the mean field functional RMF 

As already discussed and shown in Fig. ([1]) and Fig. ab-initio calculations like relativistic or non-relativistic 
Brueckner calculations have a completely different low density behavior compared to standard mean field prediction. 
For neutron matter this is a known fact and has e.g. also been investigated in the - hypothetical - unitary limit 
akp — > oo where a is the nn scattering length. In this limit many-body calculations (BHF, variational and GFMC) 
lead to a different low density behavior than RMF theory, see discussion in Ref. [23| . To account for the low density 
behavior of the DBHF equation of state one has therefore to go beyond the standard prediction of mean field theory. 
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In this paper we choose a pragmatic approach and propose a fit of the difference between the DBHF equation of 
state and the low density RMF functional. Thus we add two new functions g^ips) and g^ips) to the energy density 
so that 



eDBHF(/3«,Pp) = eRMF(Pn,Pp) + (1 - /3 )5 (pb) + P 9 (pb) , 



(18) 



where f3 = {pn ~ Pp)/pB- The additional terms g^{pB) and g^ (pb) are respectively adjusted in symmetric nuclear 
matter and pure neutron matter. The isospin degree of freedom is factorized with a quadratic function which respects 
the nuclear isospin symmetry. This approximation is often performed (see for instance but in our case, it is also 
justified afterwards by comparing the new functionals RMF+Cl and RMF+C2 to the DBHF binding energies at low 
densities (Fig. bottom panels). The functional correction is unknown, but it is clear that this correction should 
be small around saturation density, and should converge to zero at very small densities. Then an overall exponential 
shape impose to fulfill the first condition and a factorization in power of the density ensure that the second condition 
is also satisfied. We obtained two different density functionals which reproduce the data with an equal accuracy. The 
first phenomenological correction CI is a product of a polynomial function in the baryonic density with an exponential 
and has the following form: 



g'iPB) 
g'^iPB) 



Pb 
0.06 



vfPB 



(19) 
(20) 



with the parameters: t>f =-48.834 MeV.fm^, i;f =1073.3 MeV.fm^, w|=-14813 MeV.fm^, =0.03114 fm" 



,,JV_ 



MeV.fm^, p^=0.0937 fm" 



, Vq =5.373 



The second phenomenological correction C2 is a sum of two exponentials of the form 



0.01 



0.06 



g 



(22) 



P'o^ 



:0.0140 fm- 



Pi 



f =0.0879 fm- 



-0.334 



with the parameters: ■i;Q=-9.28 MeV.fm^, ■i;f=-5.48 MeV.fm^, 

MeV.fm^, =4.818 MeV.fm^, p^=0.0629 fm'^, =0.1046 fm'^. We expect to have a measure of the error induced 
by the peculiar choice of the functional by comparing the prediction obtained with the two functionals RMF-I-Cl and 
RMF-t-C2. 



Symmetric matter 



pure neutron matter 




P [fm" 



P [fm' 



FIG. 5: (Color online) Difference between the DBHF calculation and the low density RMF fit (square symbols). The corrections 
(solid lines) are drawn for symmetric nuclear matter (left panel) and pure neutron matter (right panel). 



In fig. [5l we shown the result of the adjustment of the functional CI (solid line) and C2 (dashed line) to the 
difference between the DBHF binding energy and the low density RMF functional (square symbols) in symmetric 
nuclear matter (left panel) and pure neutron matter (right matter). Despite the different density dependence between 
the functionals CI and C2, those two functionals reproduce the square symbols at the same level of accuracy. 
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D. Properties of the functionals 



The table ini give the properties of the functionals RMF, RMF+Cl and RMF+C2 around saturation density: the 
binding energy Bq, the saturation density po, the incompressibility Kq and the symmetry energies a] and a^. The 
properties of the DBHF calculation are also indicated. The properties of the low density RMF functional differ 
significantly from the DBHF results. Indeed, the saturation density and the compression modulus are higher than 
the DBHF results. It is a consequence of the low value of the parameters /^', as expected. This parameter could 
not be changed, as it is already adjusted to the quadratic density dependence of the scalar mass. However, the 
saturation properties of the corrected functionals RMF+Cl and RMF+C2 are very close to the DBHF calculation. 
Then, in our framework, the reduction of hte incompressibility modulus is induced by the low density physics. This 
is a very different understanding from the standard one which relics on non-linear corrections at high densities. We 
have also calculated the symmetry energy in two different manner: either assuming a quadratic dependence in the 
asymmetry parameter /? (aj), or as we have a functional, by performing the second derivative around symmetric 
nuclear matter (a^). The latest calculation is the exact one. We note a systematic underestimation by about 1-2 MeV 
of the symmetry energy assuming a quadratic behavior. This is a small error with respect to the difference in energy 
between symmetric nuclear matter and pure neutron matter. 



fits 


Bo [MeV] 


Polfm 


Ko[MeV] 


ai[MeV] 


ai[MeV] 


RMF 


-16.08 


0.1933 


365 


35.8 


37.7 


RMF-HCl 


-16.27 


0.1857 


251 


35.1 


36.9 


RMF+C2 


-16.24 


0.1856 


242 


35.1 


36.9 


DBHF 


-16.15 


0.1814 


230 


34.4 





TABLE II: Properties of the functionals. The symmetry energy has been obtain in two different ways: assuming a quadratic 
dependence of the equation of state between symmetric and pure neutron matter (al), or performing numerical derivation of 
the binding energy around symmetric nuclear matter (a^). A systematic difference is observed but is less than 5%. 



In this context it is worth to mention that RMF fits to finite nuclei require relatively high compression moduli 
K 300 MeV [H, [s^]. Equations of state with a stiff high density behavior stand, however, in contrast to the 
information extracted from heavy ion reactions 0,143]. The pure RMF fits to the DBHF EOS, i.e. discarding the 
low density correction term, provide equations of state which are stiff, however, not due to their high density behavior 
but do to the low density part. The compression moduli of the pure RMF contributions without correction term are 
K=365 MeV in contrast to the soft DBHF EOS with K=230 MeV. If one assumes that the correction terms contain 
effectively contributions from the deuteron and/or reflect the precursor of a superfluid low density state which leads to 
additional binding in infinite matter, but plays no substantial role in finite nuclei, this could explain the discrepancy 
between the EOS obtained from RMF fits to finite nuclei and that predicted by DBHF or the a low density virial 
expansion [3]. 



IV. SPINODAL INSTABILITIES 



Let us consider asymmetric nuclear matter characterized by a proton and a neutron densities pi = Pp, Pn- In infinite 
matter, the extensively of the free energy implies that it can be reduced to a free energy density : i^(T, V, Ni) = 
VT{T,pi). The system is stable against separation into two phases if the free energy of a single phase is lower than 
the free energy in all two-phases configurations. This stability criterion implies that the free energy density is a convex 
function of the densities pi. A local necessary condition is the positivity of the curvature matrix : 









dpidpj 




9pj . 



(23) 



where we have introduced the chemical potentials fij = ■^^\T,v,Ni — §^\T,pi^j- In the considered two-fluids system, 
the [J^ij] is a 2 * 2 symmetric matrix, so it has 2 real eigenvalues : 



± VTr[-^d -4Det 



(24) 



associated to eigenvectors Sp deflned by {i ^ j) 

6pf __ 



T 

A± - J^, 



33 



A^ — J~n 



(25) 
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Eigenvectors associated with negative eigenvalue indicate the direction of the instabihty. It defines a local order 
parameter since it is the direction along which the phase separation occurs. The eigen values A define sound velocities, 
c, by = i^Pi A. In the spinodal area, the eigen value A is negative, so the sound velocity, c, is purely imaginary 
and the instability time t is given by r = d/\c\ where d is a typical size of the density fluctuation. 

The requirement that the local curvature is positive is equivalent to the requirement that both the trace (Tr[^ij] = 
A+ + A^) and the determinant (Dct[jFy] = A+A^) arc positive 

Tt[T,j] > 0, and Det[T^j] > (26) 

The use of the trace and the determinant which are two basis-independent characteristics of the curvature matrix 
clearly stresses the fact that the stability analysis should be independent of the arbitrary choice of the thermodynamical 
quantities used to label the state e.g. {pp, p„) or (pi, ps). If Eq. (|26|) is violated the system is in the unstable region 
of a phase transition. 




FIG. 6: (Color online) Spinodal contour for the low density RMF functional, RMF, RMF+Cl and RMF+C2. 

We represent in Fig. ^ the spinodal contour for the low density RMF functional, RMF, RMF+Cl and RMF+C2. 
The contour for the RMF functional is very similar to the one obtained previously with non-relativistic interactions f^, 
H^, |4^. Meanwhile, spinodal densities are a little bit larger. Indeed, in mean-field models, it is well known that the 
spinodal density in symmetric nuclear matter scales with the saturation density [HOl . The spinodal densities we obtain 
is just due to a scaling effect induced by the saturation density which is slightly to large for the functional RMF. The 
corrections induce important modifications of the spinodal contour, especially in asymmetric nuclear matter. For an 
asymmetry parameter of about y ~ 0.4, the spinodal density, ps, is reduced by about 15-20%. 

In usual mean field calculations, it has been found that the direction of the unstable mode is still isoscalar in 
asymmetric nuclear matter As the spinodal contour calculated with the functionals RMF-I-Cl and RMF+C2 
differs significantly from the functional RMF, one could wonder what would be the consequences for the unstable 
mode: still isoscalar or possibly isovector ? We have represented our results on Fig. ([T]) where is shown the evolution 
of the unstable mode Sp„/Spp with the density for asymmetries ranging from y^O.l to 0.5. We represent only neutron 
rich matter, proton rich is easily deduced from the isospin symmetry property. An isoscalar mode is defined by 
Spn/Spp=l while a mode along y=cst satisfies Spn/Spp = y/{l — y). For the convenience of understanding Fig. ([71), 
the value ?//(l — y) is written into parenthesis in the legend of each curve. For the three functionals (RMF, RMF+Cl 
and RMF+C2) the unstable mode is included in between the isoscalar direction and the direction y=cst. However, 
the results obtained with the functionals RMF+Cl and RMF+C2 are very similar and differ from the one obtained 
with the functional RMF. Indeed, for the functionals RMF+Cl and RMF+C2, the unstable mode is less isoscalar 
than the one calculated with the functional RMF, hence, the fractionation mechanism should be less pronounced than 
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the one predicted with mean-field models (y, |43, |49| . The gas phase should then be less asymmetric than what was 
previously predicted based on mean-field calculations. 



c 
Q. 

to 



1 

0,9 
0,8 
0,7 
0,6 



q!"0,5 
CO 

0,4 
0,3 
0,2 
0,1 



1 — I — I — I — I — I — I — I — I — I — r 



- /' 



y=0.5 (1) 
y=0.4 (0.66) 
y=0.3 (0.43) 
y=0.2 (0.25) 
y=0.1 (0.11) 
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FIG. 7: (Color online) Evolution of the unstable mode 5p„/Spp with the density for asymmetries ranging from y—0.1 to 0.5. 



V. CONCLUSIONS 



We have obtained two functionals which give a very good description of the DBHF calculations in asymmetric 
nuclear matter. Those functionals are based on the low density RMF Lagrangian which is developed as a series 
expansion of kp/M. Furthermore effective terms are added which account for effects beyond mean field like the 
deuteron pole. This parametrisation has been used to understand the topological properties of the energy density, 
like the spinodal zone. We have observed that the spinodal zone is reduced in asymmetric nuclear matter contrarily 
to all the previous mean-field calculation. Based on the analysis of the direction of the unstable mode, it is shown 
that the gas phase could be less asymmetric than what was previously predicted based on mean-field calculations. 

This calculation has been performed at zero temperature while experiments probe the liquid-gas phase transition 
near the critical temperature [soj . An extension of this work to finite temperature is then necessary, but one could 
expect from mean field calculations that the critical density scales with the spinodal density ^ 0.5ps. In a future 
work, it would then be interesting to evaluate the effects of the two-body correlations at finite temperature. 

As a final conclusion and outlook of this work, we would like to stress that the low density EoS is required for 
both heavy ion collisions as well as for the description of the crust of neutron stars where a low density neutron gas 
is formed. The density functional theory is then an interesting framework where these two different nuclear systems 
could be describe by a unique functional. 
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APPENDIX A: RELATIVISTIC MEAN FIELD MODEL 

A set of coupled field equations for the meson and nucleon fields can be obtained from the Lagrangian in Eq. ([T]). 
This field equations in a mean field approximation (MFA) are 

(17^,5^ - {M - g^a - gsT^Sz) - Qu^i^ujq - gp-f"T3po)^p = 0, (Al) 

TO^cr + atJ^ + ba^ — gc < V'V' >= gaPs, (A2) 

mlivo = g^ < ip-f'^ip >= gu;PB, (A3) 

nT-lpo = 9p < -ii^l^Tsip gpP3, (A4) 

iTT-ISs ^ gs < ■4'T3ip >= gspss, (A5) 

where p^ = pp - p„ and ps3 = Psp ~ Psn, Pb ^ Pp + Pn and 

^ ' d^k-^ (A6) 



(2^)3 7o E*ik) 

are the baryon and the scalar densities, respectively. 

Neglecting the derivatives of mesons fields, the energy-momentum tensor in M FA is given by 

Tpi, = i'ipjpdi.ilj + [^TO^cr^ + U{(t) + \^m\P - ^mluixuj^ ~ ]^Tnn?pP\P^]gpu ■ (A7) 

The equation of state {EOS) for nuclear matter at T=0 is straightforwardly obtained from the energy-momentum 
tensor. The energy density has the form 

. =< T«" >= E 2 / (Strife) + + U{a) + \mlu:l + +\mlSl . (A8) 

i—7i,p 

The pressure is given by 

i—l i—n,p ^ ^ I \ y 



where E*{k) = \l k'^ + Mi*'^ , i ~ p,n. The nucleon Dirac masses are, respectively 



Mp* = M - g^cj - gsSs , (AlO) 
Mn* = M - g^a + gsSs . (All) 



In mean field approximation the kinetic contributions to energy density and pressure in Eqs. (jASp and (jA9p can easily 
be evaluated by partial integration which yields 



= E 2/ (03^*(fc) - E (Jp.i^r(fcFj + Jm;p,,) (A12) 



t—n,p 

-El/ - E (Ip^mk..) - Im:ps.) . (A13) 



(27r)3i?*(fc) 

The nucleon chemical potentials pi are given in terms of the vector meson mean fields 



pi ^ ykp- + Mi* + gujLUo T dpPo (+ proton, - neutron) , (A14) 
where the proton/neutron Fermi momenta kp. are related to the corresponding densities by fc^^ = (Stt^p^ 



U/3 
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APPENDIX B: SOLUTION METHOD OF THE NON-LINEAR SELF-CONSISTENT EQUATION FOR 

THE ISOSCALAR SCALAR FIELD a. 

Our starting point is the self-consistent equation in Let's call ctq the solution of this linear self-consistent 
equation, gaC^o = f^Psi and (Ti the first order correction induced by the non- linear terms. Then cti fulfils the following 
equation 

5(tO-i = %-(.g(TCTo + gacri f - 7 ^-^(ffo-CTo + g<jO-i)^ . (Bl) 

We first suppose that (Ti/cto ^ 1 (we will verify this hypothesis after- while). Then, Eq. (jBip leads to 

^ = ~^g.a, f 1 + 2^) - ^ {g^a,f f 1 + 3^) + . (B2) 

do mi \ ax) gam% \ ax) ctq 

Indeed, it reads 2a^x -I- ^^x^ where x — gaaQjm^ — g^ ^'m~^ P^' With typical values e.g. from the NL3 model, 
i.e. g„ ~ 10, mg- ~ 500 MeV, the parameter x is approximately x — 1.6 at saturation density. Then, the two terms 
of the denominator are about 20. The correction ai/a^ can be expressed as a function of the scalar density as 

(^0^ bat 9 

^= +o(^') (B3) 



and taking into account only first order corrections, we arrive at 

2a^ b 

gaml"^""""' ' \glm'^ {mag^) 



= -—IjI^9<y(Jof + - f_ ^ ^2 ) {ga(^of + o{g„al) . (B4) 



We keep only the first term which contributes up to the power 9 in kpi. Hence, the approximate solution of the 
non-linear self-consistent equation is 

5.^ = Os-^pl + oik%^). (B5) 
gey mi 

The first term on the r.h.s. of Eq. (jBSp is the solution of the linear self-consistent equation, while the second term is 
induced by non-linear corrections. 
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